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Summary 


A number of wire mesh and metal felt test samples, with a range of porosities, 
yield generic correlations for friction factor, Nusselt number, enhanced axial 
conduction ratio and overall heat flux ratio. This information is directed pri- 
marily toward stirling-cycle regenerator modelers, but will be of use to anyone 
seeking to better model fluid flow through these porous materials. 

Behind these results lies an oscillating-flow test rig, which measures pumping 
dissipation and thermal energy transport in sample matrices, and several stages 
of data-reduction software, which correlate instantaneous values for the above 
dimensionless groups. Within the software, a theoretical model reduces instan- 
taneous quantities from cycle- averaged measurables using standard parameter 
estimation techniques. 



Chapter 1 

Introduction 


1.1 The Continuing Saga 

This report is actually the fourth in a series of reports on progress milestones 
for this regenerator test rig. A 1989 report [3] covered test-rig design and data- 
reduction methodology, but no actual testing, under a NASA Lewis SBIR (Small 
Business Innovative Research) contract to Sun power Incorporated. The hard- 
ware has been considerably revised since then but the data-reduction methodol- 
ogy remains about the same. A 1992 report [4] presented preliminary pressure- 
drop and heat-transfer correlations for a few woven- screen and metal-felt ma- 
trices under the auspices of a NASA grant to the Ohio University Center for 
Stirling Technology Research (CSTR). And a 1994 report [13] summarized a 
larger number of tests, run but not yet completely reduced, as the conclusion 
to the test program under the now defunct CSTR. 

In this report, all of the previously-reported screen and felt experimental 
data have been reduced into generic correlations intended for application to 
all woven-screen matrices or all metal-felt matrices. This final step was made 
possible through NASA service-contract funding to Gedeon Associates. We have 
also made some effort to make this report reasonably self-contained by including 
a description of the hardware, theory and software in several appendices. One 
exception to this is the matter of error analysis, the theory of which seemed 
too lengthy and convoluted to reproduce here. While this report mentions 
various computed error estimates, it does not deal much with the theory behind 
them. Interested readers will have to refer to the original 1989 report [3] for the 
theoretical details. 


1.2 Samples Tested 

We tested a total of seven matrices, three woven-screen type and four metal-felt 
type. All samples are representative of regenerator matrices used in Stirling 
machines: 
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200 mesh 200 mesh (per inch) stainless steel woven wire screens, sintered, wire 
diameter 53.3 microns (0.0021 in), porosity 0.6232, sample thickness 10.1 
mm 

100 mesh 100 mesh (per inch) stainless steel woven wire screens, stacked, wire 
diameter 55.9 microns (0.0022 in), porosity 0.7810, sample thickness 17.5 

mm 

80 mesh 80 mesh (per inch) stainless steel woven wire screens, stacked, wire 
diameter 94.0 microns (0.0037 in), porosity 0.7102, sample thickness 22.1 
mm 

2 mil felt inconel metal felt, round wire 50.8 micron (0.002 in) diameter, poros- 
ity 0.688, sample thickness 7.54 mm 

1.5 mil felt (top sample) stainless steel metal felt, round wire 38.1 micron 
(0.0015 in) diameter, porosity 0.730, sample thickness 7.67 mm 

1.5 mil felt (middle sample) same material, porosity 0.748, sample thick- 
ness 7.49 mm 

1 mil felt stainless steel metal felt, round wire 25.4 micron (0.001 in) diameter, 
porosity 0.8200, sample thickness 12.7 mm 

0.5 mil felt stainless steel metal felt, round wire 12.7 micron (0.0005 in) di- 
ameter, porosity 0.8405, sample thickness 14.9 mm 

The metal-felt matrices were all manufactured by the Technetics (metal 
fiber) division of the Brunswick Corporation, which has since been re-organized 
into Memtec America Corporation, Fluid Dynamics Division, Deland FL, (904) 
822-8000. Memtec continues to produce fine-gauge drawn wires of round (or 
nearly round) cross section and sintered metal felts made from these wires. 
These felts have a random fiber orientation with the fibers lying predominantly 
transverse to the flow direction. The samples labeled top and middle were re- 
spectively cut from the top and middle regions of a single regenerator matnx in 
order to evaluate any difference due to manufacturing process. Evidently there 
is a small difference in porosity, probably due to the portion near the middle 
being somewhat shielded from the metal rollers used to compress the sample to 
its nominal thickness. A micrograph of a typical metal-felt matnx appears in 

figure 1.1. ... 

The screen matrices were all assembled in-house by stacking together indi- 
vidual layers in random orientation. Only one of the above samples was sintered 
because the previous report [4] concluded that there was no significant difference 
in heat-transfer or pressure-drop characteristics produced by sintering compared 
to cold stacking. In either case though, there is some variability in packing den- 
sity for the final matrix assembly, depending on the pressure exerted to compress 
the layers. Our procedure involved pressing layers together between lathe cen- 
ters to achieve a desired porosity. Quoted porosity is overall-matrix porosity 
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rather than ideal single-screen porosity. A source for woven screens is Cleve- 
land Wire Cloth, Cleveland OH, (216) 341-1832. A micrograph of a typical 
woven-screen matrix appears in figure 1.2. 
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Figure 1.2: Micrograph of a woven-screen matrix magnified lOOx. 


Chapter 2 

Screen and Felt 
Correlations 


Presented here are correlating expressions, derived from test data, for regener- 
ator pressure-drop and heat transfer in terms of mean flow variables available 
in simulation models or from experimental measurement. Symbols are defined 
in appendix A. Details of testing and data reduction are the topic of the next 
section. 

Although all samples were metallic (usually stainless steel) there is nothing 
preventing these correlations from applying to non-metallic matrices of similar 
structure. Polyester- fiber felts, in particular, are sometimes used as stirling- 
cooler regenerator matrices and are likely to be structurally identical to metal 
felts. The validity of this assertion rests on the fact that flow friction, heat 
transfer and enhanced axial conduction are all phenomena that occur solely 
within the gas as it flows around the solid obstructions of the matrix. Solid 
conduction and specific heat are independent phenomena whose affect on our 
final correlations was judged to be negligible over the experimental conditions 
encountered. Even when solid phenomena are significant (below 10 K, for exam- 
ple, when heat capacity drops off), they are properly modeled within the solid 
itself and not within the gas. Physically, the matrix solid can only affect the 
surface temperature seen by the gas, and the gas sees that surface temperature 
just like any other. 

2.1 Darcy Friction Factors 

We started out modeling friction factors in the standard two-parameter Ergun 
form [9] 

/ — ( 2 * 1 ) 

where a\ and a % were to be determined. However, we were able to get a much 
better fit to data by introducing a relatively minor modification to the Ergun 
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equation. We call this the three-parameter modified-Ergun form: 

/ = a/Re + a 2 K 3 < 2 - 2 ) 

The idea is that parameter 03 will be negative but small, allowing the correlation 
to better track observed reality at high R e . At low R e (creeping viscous flow) 
the ai/Re will dominate, exactly as in the original Ergun equation. 

Estimated parameters 01 , o 2 and 03 appear in the following tables for the 
various matrices tested: 


Woven Screens 

200 mesh 

100 mesh 

80 mesh 

a x 

02 

03 

Rem range 
V a range 
6/L range 

129.3 ±0.3 
2.990 ± 0.038 
-0.0758 ±0.0019 
0.45 — 1700 
0.0052 — 3.0 
0.028 — 2.2 

138.9 ±0.5 
2.567 ± 0.020 
-0.0816 ± 0.0010 
0.88 — 5200 
0.026 — 17 
0.087 — 1.2 

120.1 ± 0.5 

2.369 ±0.020 
-0.0836 ± 0.0010 
2.2 — 6100 
0.042 — 21 
0.079 — 0.95 

± values are 

90% confidence intervals (see appendix D) 


Metal Felts 

2.0 mil 

1.5 mil (top) 

1.5 mil (mid) 

01 

02 
03 

Rem range 
V a range 
6/L range 

120.7 ±0.5 
3.730 ± 0.033 
-0.0627 ± 0.0013 
0.83 — 2500 
0.015 — 5.2 
0.19 — 2.6 
1 mil 

214.1 ±0.6 
6.629 ± 0.068 
-0.0703 ± 0.0016 
0.67 — 1500 
0.012 — 3.7 
0.19 — 2.4 
0.5 mil 

239.1 ±0.7 
8.295 ±0.071 
-0.0903 ±0.0013 
0.36 — 1800 
0.005 — 5.6 
0.15 — 2.6 


ai 

213.8 ±0.5 

211.6 ±0.7 

a>2 

4.514 ±0.041 

4.217 ± 0.064 


-0.0705 ± 0.0014 

-0.0651 ±0.0024 

Rem r an g e 

0.21 — 1800 

0.11 —900 

V a range 

0.0056 — 5.3 

0.0021 — 1.6 

6/L range 

0.043 — 1.3 

0.102 — 2.2 

± values sure 90% confidence intervals 


We had hoped to combine the above results into two master correlations — 
one for all screens and one for all felts — by introducing a porosity dependence 
into the modified Ergun equation. The correlating expression was to be the 
four- parameter form 

f = (a l /R e + a 2 R a >)l3 a * (2-3) 


where 0 is porosity. Unfortunately this expression did not seem to fit the 
data very well. As judged by plots of modeling residuals, the fit was only 
marginally better than the above three-parameter expression, applied to the 
combined screens or combined felt data. The reason must have had to do with 
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variations in the individual samples not correlated with porosity — structural 
dissimilarities at a microscopic level, for example. Even without explicit poros- 
ity dependence, porosity does affect the calculation of hydraulic diameter dh 
and void-average flow velocity u , in terms of which friction factor is defined. 
This seems to be as much porosity dependence as our data justifies. 

So, the best we were able to do was apply the three-parameter friction factor 
to the combined data sets of all screens and all felts. In doing so we lost a good 
deal of our precision of fit to the data. Reduced individually, worst-case residuals 
(deviations of data from model as a fraction of estimated random error) were 
on the order of ± 5. Reduced together, worst-case residuals were on the order 
of ± 30 for screens and ± 80 for felts ( see figure 3.1). Were it not for the 
extremely high signal to noise ratio of about 300 (see figure 3.2), this would be 
bad news indeed. As it is, it suggests a worst-case relative pressure drop error of 
about 27% for felts and 10% for screens when applying the following combined 
correlations to a generic woven-screen or metal- felt matrix (see section 3.1): 

Woven Screens 

/ = 129/Re + 2.91R- 0 103 (2.4) 

Metal Felts 

/ = 192/Re + 4.53R” 0 * 067 (2.5) 


Or in tabular form: 



Woven Screens 

Metal Felts 

ai 

129.3 ±0.2 

192.4 ±0.2 

<*2 

2.913 ±0.013 

4.533 ±0.020 

03 

-0.1027 ± 0.0006 

-0.0672 ± 0.0007 

R,_ m range 

0.45 — 6100 

0.11 —2500 

V a range 

0.0052 — 21 

0.0021 — 5.6 

8/ L range 

0.028 — 2.2 

0.043 — 2.6 


± values are 90% confidence intervals 


The above correlations are plotted in figures 2.1 and 2.2. 
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Re 

Figure 2.1: Oscillating-flow friction factor for the combined woven screen test 
samples 
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2.2 N u and N k Simultaneously 

Nusselt numbers and enhanced conductivity ratios as presented here should 
be used together in an oscillating-flow regenerator model that simultaneously 
includes the effects of both film heat transfer and enhanced axial conduction. 
Since they are based on best estimates of the true regenerator physics, they 
may also be used with caution by porous flow modelers in other fields. See 
appendices C and F for further discussions on the nature and use of N u and N k . 

Results for the individual woven-screen and metal-felt samples are not pre- 
sented since, this time, we were more successful in combining data sets to pro- 
duce two master correlations — one for all screens and one for all felts. The 
likely reason for our success is that porosity has a fundamental physical affect 
on Nusselt number N u and axial conductivity ratio N k that overwhelms any 

structural variations. . 

The presumed correlating expressions are in terms of Peclet number ( RtP r ) 

and porosity: 

N u = (1 + a l P?W 3 (2-6) 

N k -N k o = a 4 i?’0 O8 (2 ' 7) 

The constant 1 in the Nusselt number expression is a concession to prevent 
N u from vanishing for P t = 0, which may be important in preventing ther- 
mal decoupling of the gas and matrix, and possible instability, in numerical 
simulations. It is small enough not to significantly affect our parameter estima- 
tion. Some constant value is justified by laminar solutions in simple geometries, 
where N u tends to some fixed non-zero value N u0 as P t approaches zero. A 
value of N u0 « 1 seems reasonable in light of the published limiting value of 
hd w /k = 0.43 (Nusselt number based on wire diameter) in [6] for flow normal to 
a single cylinder. The definition of N u employed in the present report is based 
on hydraulic diameter and is a factor of 0/(1 ~ 0 ) ^ger. Unfortunately, due to 
the dominance of enhanced conduction as the axial energy transport mechanism 
for low P e 's, it is not possible to estimate N u0 directly from our data. 

That the P e exponents for both N u and N k are the same reflects the as- 
sumption that the underlying enhancement mechanisms are the same for both. 
In essence, the argument is that fluctuating parts of velocity and temperature 
fields enhance axial thermal diffusion by the same factor they enhance gas-to- 
matrix heat transfer. This topic is picked up again in appendix F. It would be 
nicer if we didn’t have to make this assumption but previous efforts have shown 
that it is impossible to independently correlate P e dependence for both N u and 
N k from our data. 

On the other hand, it does turn out to be possible to independently correlate 
porosity dependence for N u and N k . Based on considerations presented below, 
it is arguable that the asymptotic form for Nusselt number as 0 -> 1 should scale 
as — /?), raised to some positive exponent while the asymptotic form for 
enhanced axial conductivity should do likewise, except for a negative exponent. 
But when we tried correlating N u and N k in terms of (0/(1 - 0))“ or even 
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(1 — /3)°, (exponents a different for N u and N k ) we found the fit to the data, 
as measured by the residuals, was significantly worse than when correlated in 
terms of /?°. Evidently, our porosities were never high enough for the asymptotic 
behavior to kick in? Highest porosities tested were 0.78 for screens and 0.84 for 
felts. So, on empirical grounds we must say that over the limited range of 
porosities represented in our samples, a porosity dependence in the form (3 0, 
appears to fit the data best. 

Final correlating expressions for the two types of material are: 


Woven Screens 


N u 

N k - N k0 

Metal Felts 

N u 

N k - N k o 


(1 + 0.99 if • 66 )^- 7 » (2.8) 

0.50 ffiKp- 2 - 91 (2.9) 

(1 + 1.1 e ^? 66 )/? 2 - 61 (2.10) 

1 30 ^. 660 - 2.°9 ( 2 . 11 ) 


Or in tabular form: 



Woven Screens 

Metal Felts 

ai 

0.991 ± 0.048 

1.159 ±0.033 

a 2 

0.662 ± 0.005 

0.656 ±0.005 

<*3 

1.792 ±0.048 

2.609 ±0.049 

o 4 

0.504 ±0.12 

1.299 ±0.090 

05 

-2.908 ±0.48 

-2.089 ±0.24 

Rtm range 

1.04 — 3400 

0.79 — 1400 

V a range 

0.0048 — 16 

0.0037 — 3.3 

<5/L range 

0.17 — 3.0 

0.17 — 3.8 


± values are 90% confidence intervals 


The above correlations are plotted in figures 2.3 and 2.4 and residuals are 
plotted in figure 3.3. 

One cannot help but notice the strong dependence of both N u and N k on 
porosity and that N u increases with porosity while N k decreases. That the 
dependence is strong is tempered by the observation that porosity typically 
only varies over the range of about 0.60 to 0.90 in most regenerators. Even the 
strongest porosity factors in the above correlations vary only by a factor of about 
3 over this range. Perhaps more interesting are the trends of the correlations. 
These are at least in the correct direction, which we can deduce by arguing 
about the high-porosity limiting case where the matrix looks like empty space 
with a few widely scattered wires. For then, Nusselt number (based on hydraulic 
diameter) will tend to the value for a single cylinder in cross flow, which has 
already been remarked to scale by the ratio of hydraulic diameter to wire diam- 
eter, or (3/(1 - /?), which increases as /? increases. This is somewhat countered 
by the Reynolds number dependence of single-cylinder Nusselt number which 
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scales roughly as Re° w 5 over the range of interest, where R ew is Reynolds number 
based on wire diameter, rather than hydraulic diameter. When cast in terms 
of our Peclet number this introduces a factor (0/(1 - /?))“ 0 ' 5 into N u , enough 
to reduce but not overwhelm the previous factor. In the case of enhanced axial 
conductivity, we expect Re w to again characterize the flow in the limit /? — ► 1 
so that Nk should also scale roughly as Rj ^ 5 . We might also throw in a fac- 
tor of (1 - /?)//?, the filled- to- void volume ratio, to account for the fraction of 
the volume filled by turbulent wakes. Based on these considerations, we expect 
a porosity dependence of roughly ((1 - £)//?) 15 in our N k expression, which 
decreases as /3 increases. 
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Figure 2.3: Oscillating- flow Nusselt number and enhanced axial conductivity 
ratio for woven screen, for (3 = 0.75 





Nko @ /? = 0.75 N u @ (3 = 0.75 
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2.3 N ue assuming Nk = Nko 

Effective Nusselt number N ue captures the total oscillating- flow regenerator 
heat flux ql in enthalpy transport alone. We intended it for use primarily by 
Stirling modelers who prefer to neglect enhanced axial conduction. In other 
applications it will tend to under-predict actual film heat transfer for a given 
temperature difference — especially at Peclet numbers below about 10. The 
presumed correlating expressions are 

N ue = (1 + a 1 P < ?*)/?* 3 (2.12) 

N k -N k o = 0 (2.13) 

The constant 1 in the Nusselt number expression and the form of the porosity 
dependence are there for the same reasons as previously discussed for simulta- 
neous N U1 N k modeling. 

Final correlating expressions for the two types of material are: 


Woven Screens 

N ue = (1 + 0.64 P? n )0 179 

(2.14) 

Metal Felts 

N ue = (1 + 0.48P®' 79 )/? 2 ' 75 

(2.15) 


Or in tabular form: 

Woven Screens Metal Felts 


ai 

0.644 ±0.017 

0.485 ±0.007 

d2 

0.720 ±0.003 

0.794 ±0.002 


1.794 ±0.030 

2.752 ±0.028 

Rem range 

1.04 — 3400 

0.79 — 1400 

V a range 

0.0048 — 16 

0.0037 — 3.3 

6/L range 

0.17 — 3.0 

0.17 — 3.8 


± values are 90% confidence intervals 


The above correlations are plotted in figures 2.5 and 2.6 and residuals are 
plotted in figure 3.4. The narrower confidence bands this time reflect the com- 
parative ease of modeling three rather than five parameters, not any intrinsic 
superiority of the world view neglecting enhanced axial conduction. In fact, the 
previous five-parameter modeling for simultaneous N u and N k gave a slightly 
but significantly better fit to the data, as measured by the minimum chi-squared 
value of the residuals. 



CHAPTER 2. SCREEN AND FELT CORRELATIONS 


17 



Figure 2.5: Effective Nusselt numbers neglecting enhanced axial conduction, for 
woven screen, for (3 = 0.75 
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Figure 2.6: Effective Nusselt numbers neglecting enhanced axial conduction, for 
metal felt, for P = 0.75 
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2.4 Overall Heat Flux Ratios 


Overall heat flux ratio correlates the overall thermal loss down an oscillatmg- 
flow regenerator with a simple expression in terms of the peak Peclet number 
and porosity. 

N q = a x P € %P* 

Note that Peclet number is the peak not instantaneous value. As judged by the 
minimum chi- squared value of the residuals (figure 3.5), this expression does 
not fit the data as well as the previous expression for N u and N k , but somewhat 
better than the expressions for N ue . This is a rather remarkable achievement 
for such a simple idea. 

Final correlating expressions for the two types of material are: 


Woven Screens 
Metal Felts 


N q = 0.194P e ^ 30 /? -1 81 

Nq = 0.253Peii 24 r 267 


(2.16) 

(2.17) 


Or in tabular form: 



Woven Screens 

Metal Felts 

ai 

0.194 ±0.005 

0.253 ± 0.003 


1.301 ±0.003 

1.239 ±0.002 

a 3 

-1.806 ±0.030 

-2.672 ± 0.029 

R cm range 

1.04 — 3400 

0.79 — 1400 

V a range 

0.0048 — 16 

0.0037 — 3.3 

6/ L range 

0.17 — 3.0 

0.17 — 3.8 


± values are 90% confidence intervals 


The above correlations are plotted in figures 2.7 and 2.8 and residuals are 
plotted in figure 3.5. 





N v @ p = 0.75 




Chapter 3 

Test Details 


3.1 Pressure-Drop Tests 

Pressure-drop testing is relatively straight forward, involving only a measure- 
ment of sample pumping dissipation inferred from the PV work done by the 
piston. The principal uncertainties are due to the compressibility of the work- 
ing gas when pressure-drop is large. 


3,1,1 Data Sets 

Individual samples were tested as follows: 

Pressure Drop Tests 

test dates no. points 


200 mesh screens 

11-03-92 

169 

100 mesh screens 

11-04-92, 11-05-92 

237 

80 mesh screens 

11-02-92 

193 

2.0 mil felt: 

6-29-93 

176 

1.5 mil felt (top) 

6-30-93 

148 

1.5 mil felt (mid) 

6-02-93 

159 

1.0 mil felt 

11-04-92 

201 

0.5 mil felt 

11-03-92 

168 


Derived data for all screens were collected into file PSCN.DRV and for all felts 
into PFELT.DRV as the basis for the final correlations in this report. 

3,1.2 Error Analysis 

In data modeling parlance, residuals measure the fit of our model to the exper- 
imental data and show the nature of any discrepancy between the two. This 
information cannot be directly plotted on our friction factor curves because our 
experimental measurable is not instantaneous pressure drop but, rather, cycle- 
mean sample pumping dissipation. So we need separate residual plots, which 
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the data modeling software produces each time it runs. Plotted residuals are 
normalized by the estimated random measurement errors, which in this case 
are the random errors expected in calculating cycle-mean pumping dissipation, 
a quantity derived for each data point in terms of position and pressure trans- 
ducer errors. The data modeling software also produces plots of signal- to- noise 
ratio, which in this case show the actual values of measured cycle-mean pumping 
dissipation, normalized by the same estimated random measurement errors. 

Residual plots for our master screen and felt friction factors are shown in 
figure 3.1. Ideally, residuals should fall randomly between ± 1, when the data 
fit the model. With actual residuals falling in the range ± 30 for screens and 
± 80 for felts, this is decidedly not the case. This can mean one of two things: 
Either the model is bad or the random error is so small that even a small 
discrepancy in the model shows up as a large residual error. In the present case 
the latter phenomenon appears largely responsible. The two measurements that 
contribute to random measurement error are time-varying piston position and 
pressure. Estimated measurement errors for these were 36 microns (0.0014 in) 
and 2.8 Pa (0.0004 psi), respectively — remarkably small. These estimates were 
verified by analysis of data on a relatively small test sample. Remember that 
the errors we are talking about here are the random errors produced by two 
separate readings of the same piston position and pressure, rather than the 
systematic error, in absolute terms, due to mis-calibration or whatever. The 
true impact of these small errors may be seen in the plot of signal-to-noise 
ratios in figure 3.2 which range up to about 300 for both screens and felts. A 
residual of 80, presuming it corresponds to a high-signal data point, is therefore 
only about 80/300 = 0.27 of the total signal while a residual of 30 is only about 
0.10 of the signal. And these are worst-case errors. From this we can expect 
our correlations to predict pressure drop within about 27% for felts and 10% 
for screens, at worst. Agreement will probably be better for most cases. 

Although not shown, residuals for individual-sample friction factors are much 
lower. In fact they fall in the range of ±5, indicating a much better fit of data 

to model. . 

One additional source of error worth considering in pressure drop testing is 
the error in the model-predicted pumping dissipation as a result of variations in 
mass flux g across the test sample. Mass flux varies because of time-dependent 
density variations resulting from the pressure swing in the cylinder volume. The 
pressure-drop data modeling software produces an error plot which estimates 
the relative magnitude of this error. This error was small for all samples tested. 

3.1.3 Comparison to Other Data 

Published data for friction factors in steady porous flow is relatively common. 
Osdllating-flow results are much rarer, but this does not worry us too much be- 
cause of the negligible role Valensi number (measure of dimensionless frequency) 
seems to play in our data, as was demonstrated in (4). 

For comparison purposes we used results reported for screens in figure 7-9, 
p. 149, of Kays and London [7], Because of the lack of an analytic formulation 
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Figure 3.1: Residuals for modified- Ergun friction- factor parameter estimates: 
combined woven screens \ 2 — 7.22E4 with 599 points (top), combined metal 
felts x 2 = 8.62E5 with 852 points (bottom). 
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Figure 3.2: Signal to noise ratios for friction-factor parameter estimates: com- 
bined woven screens (top), combined metal felts (bottom). 
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for the Kays and London data, we have not plotted a comparison of our woven- 
screen friction factor against theirs. Rather, we have tabulated the agreement at 
several representative Reynolds numbers. When comparing our Darcy friction 
factors to Kays and London’s Fanning factors we must always remember to 
multiply by 4. 


Screen Friction Factor Comparison 


Re 

8 

40 

100 

400 

1000 

6000 

f — low Kays & London 

18 

5.2 

3.0 

1.8 

1.4 

1.1 

/ — present report 

18.5 

5.21 

3.10 

1.89 

1.56 

1.40 

/ — high Kays & London 

29 

8.5 

4.5 

2.7 

2.1 

1.6 


The table covers a Reynolds number range from 8 to 6000, which was the range of 
overlap for our respective data sets. Peak Reynolds number for our experiments 
ranged from 0.45 to 6100, while Kays and London steady Reynolds numbers 
ranged from 8 to 10 5 . It is seen that our woven-screen friction factor falls on 
the low side of the Kays and London data for low Reynolds number while near 
the middle for high Reynolds numbers. The Kays and London data was based 
on four discrete porosities: 0.832, 0.766, 0.725, 0.602, with either 0.832 or 0.602 
(depending on Re) producing the high friction factor and 0.725 producing the 
low friction factor. Interestingly, their data, like ours, shows no clear trend 
with porosity. And, there is about a ± 25% scatter among their four samples, 
in the ballpark of our ± 10% accuracy projections based on relative pumping- 
dissipation error. 

The Kays and London correlations were “established from tests of wire 
screens at the low Reynolds number end but crossed-rod matrices at the high 
Reynolds number end” , as they explain it. Data for screens alone was reported 
in Coppage and London [2] for Reynolds numbers up to about 1000 and the 
same porosities appearing in Kays and London [7], Apparently, the Coppage 
and London data were combined with later crossed-rod data to produce the Kays 
and London correlations. But, the structure of a crossed-rod matrix (comprised 
of straight rods) is noticeably different than a true woven screen; a fact that 
should be kept in mind when comparing our data to theirs. 


3.2 Heat-Transfer Tests 

As its principal experimental variable, heat transfer testing measures net heat 
flux down the regenerator inferred from cooler heat rejection less PV power 
and heat leaks. Compared to pressure-drop testing, this is more complicated, 
time consuming, and noisier; the noise being chiefly due to the rather large and 
time- varying nature of heat leakage. 

3.2.1 Data Sets 

Individual samples were tested as follows: 
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Heat Transfer Tests 



test dates 

no. points 

200 mesh screens 

11-17-92, 11-18-92 

273 

100 mesh screens 

11-16-92, 11-17-92 

231 

80 mesh screens 

11-12-92 

189 

2.0 mil felt: 

11-21-92, 11-22-92 

227 

1.5 mil felt (top) 

11-23-92 

272 

1.5 mil felt (mid) 

3-10-93, 3-11-93 

283 

1.0 mil felt 

11-18-92, 11-19-92 

286 

0.5 mil felt 

10-07-92, ..., 10-12-92 

271 


Derived data for all screens were collected into file HSCN.DRV and for all felts 
into HFELT.DRV as the basis for the final correlations in this report. 


3.2.2 Error Analysis 

For our heat-transfer tests we have three sets of residuals corresponding to the 
simultaneous N u and N k modeling, N ue modeling or N q modeling. In all cases 
the measured experimental variable is cycle average heat flux down the ma- 
trix, less static thermal conduction, and the normalization error is the random 
component of this, calculated for each data point in terms of various transducer 
errors. The theoretical experimental variable varies though, depending on which 
model is being compared to the data. 

The various residuals are plotted in figures 3.3, 3.4 and 3.5. Compared 
to the friction-factor residuals, these residuals are much smaller, falling in the 
range of ± 4 for both screens and felts, suggesting a better fit to the data, or 
higher noise since the signal-to-noise ratios in figure 3.6 range only to about 40. 
By comparing normalized residuals to signal-to-noise ratios at various Reynolds 
numbers, we are led to expect a prediction of overall regenerator energy flux to 
within ± 10% near the high-end of the Reynolds number range with accuracy 
falling off for lower Reynolds numbers to no better than ± 50% below Reynolds 
number about 5, using any of our correlations. Another feature of the heat- 
transfer residuals compared to the friction- factor residuals is that they are more 
randomly distributed, showing hardly any systematic deviations from the model. 
This suggests the residuals are indeed produced by random measurement errors, 
as they are supposed to be, even though these errors may be relatively large. 

The heat-transfer data modeling software also produces a number of other 
error plots, namely 

1. Relative error induced by g error 

2. Relative error induced by T error 

3. Relative error induced by neglected terms 

The first type of error is the effect on modeled regenerator heat flux of estimated 
errors in mass flux g. Both g and its error estimate come from numerical sim- 
ulation — g being the mid-sample value and its error being the variation in g 




Figure 3.3: Residuals for simultaneous N Uy Nk parameter estimates: combined 
woven screens x 2 = 1.86E3 with 693 points (top), combined metal felts x 2 — 
3.69E3 with 1339 points (bottom). 













Figure 3.5: Residuals for overall heat flux N q parameter estimates: combined 
woven screens x 2 — 2.14E3 with 693 points (top), combined metal felts x 2 = 
4.79E3 with 1339 points (bottom). 
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across the sample. The second type of error is the effect on modeled regenerator 
heat flux of estimated errors in the presumed- constant sample midpoint tem- 
perature gradient All models predict linear variation of regenerator heat 
flux with |£. The error estimate for is based on instrumentation error and 
local temperature variations related to solid- to-gas heat capacity ratio and flow 
tidal amplitude. The third type of error is the effect on modeled regenerator 
heat flux of key neglected terms in the model such as: the component of solid 
temperature variation that contributes to net enthalpy flux, time-varying pres- 
sure, spatially-varying mass flux and spatially-varying axial conduction. These 
errors only pertain to models correlating N u , Nk or N ue , not to those corre- 
lating overall heat flux N q . Errors from all these sources were uniformly small 
(combined together no more than 20% of the estimated measurement error) so 
they have no significant effect on the results and are not shown. 


3.2.3 Comparison to Other Data 


Good published data for heat-transfer properties in porous materials is scarce. 
Our only comparison is again to the data of Kays and London [7] for woven- 
screen matrices, figure 7-8, p. 148. In the low Reynolds number range (for true 
woven- screen matrices) their correlations were based on the transient- flow tech- 
nique which tracks a step temperature change in the fluid (introduced upstream 
of the matrix) as it emerges downstream under steady flow conditions. Data 
reduction correlates the heat-transfer coefficient with the maximum slope of the 
exit temperature-time plot, neglecting axial thermal diffusion. We are willing 
to ignore the differences between steady and oscillating flow because of the low 
Valensi numbers of our tests. But ignoring axial thermal diffusion is another 
thing. The presence of thermal diffusion tends to reduce the maximum slope of 
the exit temperature profile. Therefore, the maximum-slope method tends to 
underpredict heat transfer coefficient when the temperature slope is steepest — 
or, as our results suggest, for Reynolds numbers below about 10. 

We have two candidates for comparison to the Kays and London data, our 
best-guess value for the true Nusselt number N u and the effective Nusselt num- 
ber assuming zero axial thermal diffusion N ut . We shall compare both. Since, 
the assumptions behind the Kays and London correlations are consistent with 
our method of reducing N ucy it would seem that N ue offers the most meaningful 
comparison. However, both our N ue and Kays and Londons data are likely to 
suffer from error at low Reynolds (Peclet) numbers, and the errors may not be 
the same. 

In making the comparison, we should also note that Kays and London cor- 
relate heat transfer in terms of Stanton number, which is related to our Nusselt 


number by 



(3.1) 


The following tables assume P T =0.7 in converting Stanton to Nusselt number 
and Reynolds to Peclet number. 
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Screen Nusselt Number Comparison — /? — 0.602 


R e 


4 

10 40 

100 

400 

1000 

4000 

N u — 

Kays & London 

1.04 

1.81 4.16 

7.10 

16.1 

28.1 

65.6 

N u — 

present report 

1.19 

1.85 4.00 

7.00 

16.9 

30.5 

75.6 

N ue - 

- present report 

0.95 

1.45 3.25 

5.90 

15.3 

29.3 

78.7 

Screen Nusselt Number Comparison - 

0 = 

0.832 



Re 


40 

100 400 

1000 

4000 



N u — 

Kays & London 

9.72 

14.5 28.4 

46.9 

101 



N u - 

present report 

7.14 

12.5 30.1 

54.5 

135 



N ue - 

- present report 

5.79 

10.5 27.3 

52.2 

140 




The tables covers a Reynolds number range from 4 to 4000, which was the 
range of overlap for our respective data sets. Peak Reynolds number for our 
experiments ranged from 1.0 to 3400, while Kays and London steady Reynolds 
numbers ranged from about 4 (for 0 = 0.602) to 10 s (for 0 = 0.832). Although 
above 1000, as mentioned above, the Kays and London data was probably for 
crossed-rod matrices. We expect good agreement among all three Nusselt num- 
bers for high Reynolds numbers, and this is somewhat borne out by the tables. 
Best agreement is for Reynolds number about 500 with increasing disagreement 
at 4000, possibly because it falls just outside the top range of our data or pos- 
sibly because the heat-transfer characteristics of crossed rods are a bit different 
than woven screens. At low Reynolds numbers where we expect the worst agree- 
ment we see that both N ue and the Kays and London value are smaller than 
our true N u , but not by much. In fact there is remarkably good agreement 
down there, all three predicting a suspiciously low Nusselt number. This is not 
surprising in light of the likelihood that enhanced axial conduction is interfering 
with proper Nusselt number reduction for all three at low Reynolds numbers. 



Chapter 4 

Concluding Remarks 


The achievements of this effort have been both practical and theoretical. On the 
practical side we have produced correlating expressions for pressure drop and 
heat transfer for the two types of regenerator material commonly used in Stirling 
engines and coolers: woven screens and metal felts. While steady-flow data for 
woven screens has been previously reported in [7], we have confirmed its validity 
for Stirling use and extended its range to the low Reynolds numbers commonly 
found in Stirling regenerators. Data for felts is new. In a break with tradition, 
we have placed enhanced gas conductivity on an equal footing with film heat 
transfer and successfully correlated expressions for both simultaneously. And 
we have introduced the concept of overall heat flux ratio, correlated in terms of 
peak Peclet number (peak i2eP r ), which fits the data remarkably well. 

Accuracy for our derived correlations was limited by inter-sample variations 
not correlated with porosity (in the case of friction factor) and thermal noise at 
low Reynolds numbers (in the case of heat-transfer correlations). For friction 
factors, relative accuracy was a fairly uniform worst-case of 10% for screens and 
27% for felts, over the entire Reynolds number range. This should be understood 
as the worst likely relative error in cycle-average pumping dissipation that would 
be produced by our correlations for a typical sample matrix. For heat transfer 
correlations the relative accuracy was more dependent on Reynolds number, 
principally because the measured signal was always low at low Reynolds number. 
Thus the worst likely relative error in cycle-average total axial heat transport 
ranged from about 10% at peak Reynolds number on the order of 1000 to about 
50% for peak Reynolds numbers below about 5. It would be prudent then to 
avoid placing too much faith in our heat-transfer correlations for a situation 
where peak Reynolds number was much below 5. Fortunately, it is the author’s 
observation that most stirling-cycle regenerators operate in a range significantly 
above this. This is corroborated by Seume and Simon ([12], fig. 7, p. 537) who 
found no peak Reynolds number below about 30 in a survey of about a dozen 
Stirling engine regenerators in use at the time. 

On the theoretical side we have managed to dispel some of the myths sur- 
rounding oscillating flow. Under the range of conditions we tested, which were 
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intended to be representative of most Stirling applications, we found no essen- 
tial differences compared to steady flow. (The 1992 report [4] addresses this 
issue in more detail.) In other words, instantaneous local Reynolds number or 
Peclet number (/?«P r ) appear to characterize the flow quite adequately. This 
means, for example, that porous-material friction factors obtained by others un- 
der steady-flow conditions should suffice for oscillating-flow conditions in Stirling 
regenerators. The exception to this might be for Valensi numbers above about 
20, the upper range of the values we encountered. 

It is not quite so simple with heat transfer. While we agree, in principle, that 
heat-transfer and enhanced axial conduction also depend only on local instan- 
taneous Reynolds number, the trick is in reducing the data how one converts 
experimental measurables to Nusselt numbers, etc. There is no experimental 
apparatus we know of that directly measures the heat transfer coefficient in re- 
alistic regenerator matrices. The wires and void spaces are simply too small to 
introduce temperature transducers. So there is always some intermediate model 
in terms of which heat-transfer coefficient is deduced. In our case the model is a 
bit complicated as it involves total axial energy transport under oscillating flow 
conditions. Other approaches tend to track exit temperature profiles, usually 
under steady flow conditions. Both tend to neglect phenomena that compli- 
cate the analysis. The question is: which method is most likely to produce the 
correct result in a Stirling engine or cooler simulation? We would argue that 
our method does, because it is most likely to reproduce the quantity that really 
matters: total energy transport down the regenerator. 



Appendix A 

Symbols 


The following symbols are those appropriate for a computational model of 
porous regenerator flow where the computational grid is large compared to the 
matrix pore size. Accordingly, thermodynamic variables like T, u, etc., are 
understood in a local-average sense — spatially averaged over a large enough 
volume to remove eddy fluctuations, yet still small compared to the overall prob- 
lem dimensions. This viewpoint is consistent with the measurement resolution 
in the test rig itself. 
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Cp gas specific heat at constant pressure 
dh hydraulic diameter: dh, = 4 / s 

ei error of zth data point; difference between measured yt and theo- 
retical value 

g mass flow rate per unit void area: pu 
g m peak value of g 

h heat transfer coefficient between gas and solid matrix 
h e effective heat transfer coefficient that includes the effects of k a when 
used in an oscillating- flow regenerator model 
k molecular gas conductivity 

k a apparent agitated-flow gas axial conductivity (based on mean void 
flow area) 

L test-sample length 
P pressure 

qh part of Qt attributed to enthalpy transport (finite h) 

q k part of q t attributed to apparent axial conduction 

q t time-average total energy flux (per unit void area) down the 
regenerator 

q * q ty less static conduction 

s matrix surface area per unit void volume 
t time 

T gas temperature (local void average) 

T s solid matrix temperature (local solid average) 
u gas velocity (local void average) 

x axial coordinate 

y i measured value of zth data point 
ft porosity (void volume / total volume) 

/x molecular gas viscosity (local void average) 

p gas density (local void average) 

d advanced estimate of e* based on random measurement error 
(jj angular frequency (rad/s) 
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Dimensionless Groupings 

/ Darcy friction factor: dk^/(p^^) 

Nk axial conductivity enhancement ratio: ^ 

Njco value of Nk at static flow conditions 
N u Nusselt number: 

N ue effective Nusselt number: 

N q overall heat flux ratio: 

P e Peclet number: R^P r 
P em peak Peclet number: R em P r 
P r Prandtl number: 

Rc Reynolds number: ^ 

Rz m peak Reynolds number: 

V a Valensi number: 

6/L tidal amplitude ratio: Q m /{pujL) = {dh/(4L))Rc m /V a 

(3 porosity: void volume / total volume 

Subscripts 
o static-flow value 

m peak value (amplitude of first harmonic) 

Operators 

{ ) time average over one cycle: (/) = fdt 



Appendix B 

Hardware 


In its heat-transfer mode, the test rig measures exactly the quantity that is of 
bottom-line importance in Stirling machines: net thermal energy flux. And it 
does this over a wide range of operating conditions typical of Stirling designs. 

The essential idea is to attach, in sequential fashion, a cooler, regenerator 
-sample and heater to a piston cylinder. The assembly is oriented so the heater 
outlet opens into a relatively large fixed volume. When the rig is motored, an 
oscillating gas flow is set up in the regenerator sample. Moreover, due to the 
large volume of the system, the pressure swing in the piston cylinder is due 
predominantly to the frictional pressure drop of the regenerator sample and 
heat exchangers and not to compression effects. Therefore, the gas mass flow 
rate is very nearly sinusoidal and spatially uniform in all three heat exchangers. 

The rig is motored with the heater temperature elevated about 200 C above 
the cooler temperature causing a temperature gradient to develop in the regen- 
erator. As a consequence of axial conduction and imperfect heat transfer, a 
net thermal flux occurs across the regenerator, and is ultimately rejected in the 
cooler, where it is measured. Piston PV power (due to the pressure drop across 
the system) and static conduction losses through cylinder walls and flanges also 
contribute to cooler heat rejection, but these may be calibrated out. By doing 
so, the thermal performance of the regenerator in isolation can be inferred. 

Refinements of these essential ideas are necessary to reduce the above con- 
ceptual scheme to practice. Some key refinements are listed below. 

• The volume into which the heater gas flows is thermally insulated from 
its cold surroundings to prevent the influx of cold gas during the suction 
stroke from placing a large thermal load on the heater. 

• To minimize the effects of non-regenerator thermal losses, the regenerator 
sample canister is fabricated out of insulating material. Moreover, the 
cooler is also insulated from its surroundings to prevent heat from un- 
wanted sources being rejected in the coolant. The heater exterior is also 
insulated to prevent undue heating of the gas within the pressure vessel. 
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• Jet penetration from the heat exchangers into the faces of the regenerator 
might pose a problem in some samples. To avoid this as much as possible, 
flow diffusers are incorporated between the heat exchangers and the test 
sample. 

For purposes of pressure drop testing the rig is somewhat simplified. The 
cooler, heater, and flow diffusers are removed and the sample canister attached 
directly to the piston cylinder. 


B.l Early History 

The current test rig is actually a modified form of an earlier pressure-drop test 
rig designed by Sunpower Inc. for NAS A- Lewis. Reference [8] describes this 
rig in detail and figure B.l shows it in cross section. Several components of 
the earlier rig remain the same. At the heart of the rig is a variable stroke 
and variable frequency linear drive motor. A displacement section, consisting 
of a single close fitting piston in a cylinder, is directly attached to this linear 
drive motor. The test section is connected to the other end of the displacement 
section. The assembled motor, displacement, and test section are then enclosed 
in an outer pressure vessel to allow for operation at elevated mean pressures. 

The versatility of the basic test rig is shown by the wide range of operating 
parameters listed in the following table: 


stroke (cm) 
frequency (Hz) 
mean pressure (bar) 
motor power (kW) 
piston diameter (cm) 


Rig Operating Parameters 
0-3.0 
<1 - 120 

0 - 151 (absolute) 

0-2 (at 60 Hz, 3 cm stroke) 

1.9 (standard configuration, 

smaller pistons possible) 


larger and 


For regenerator heat transfer testing, the test section of the original pressure 
drop rig (a single element) was replaced by a multiple element unit consisting 
of a cooler, regenerator, heater and buffer volume. The design and fabrication 
of this setup was performed by Sunpower Inc. for NAS A- Lewis under contract 
NAS3- 25620. A detailed description of this hardware is given in reference [3]. 


B.2 Current Rig Description 

The rig was again modified during testing at the Ohio University Center for 
Stirling Technology Research (CSTR). Figure B.2 shows a cross sectional draw- 
ing of the test section of the rig as presently configured for heat-transfer testing. 
Pictures of the test rig are presented in figures B.3, B.4 and B.5. 
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Figure B.2: Cross section of current regenerator test rig in heat-transfer testing 
mode. 
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Figure B.3: Assembled Regenerator Test Rig. 
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Figure B.4: Regenerator test rig internal assembly with insulation installed. 
Scale shown by co-author Gary Wood 
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Figure B.5: Regenerator test rig internal assembly detail, without insulation. 
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B.2.1 Rig Orientation 

The rig had been originally designed for operation with the hot end down. 
Early testing at CSTR indicated that significant convection losses were occurring 
in spite of the installed insulation. So, the rig was inverted to reduce these 
convection losses. This also allowed the test rig to come to steady state more 
quickly, especially.when the pressure level in the rig was changed. 

B.2.2 Buffer Volume 

The buffer volume, attached to the end of the heater opposite the end connected 
to the test section, is a thermally-insulated canister for the purpose of containing 
the hot gas exiting the heater. An inadvertent side- effect was that it prevented 
timely pressure charging or discharging because the only gas pathway from the 
buffer volume to the pressure vessel interior was the piston clearance seal, which 
was very effective. 

So, we added the capillary tube seen in the drawing to vent the buffer volume 
to the interior of the outer pressure vessel, sized so as to give reasonable charging 
times, while not affecting the accuracy of the measurements. 

B.2.3 Heater Section 

The heater section was originally a shell and tube heat exchanger. Heat was 
supplied by means of a pumped (closed) Dowtherm loop. The heating element 
and the pump were external to the outer pressure vessel of the test rig, with the 
fluid carrying lines passing into the vessel through the pressure vessel extension. 

To simplify this system we replaced it with a copper cylinder having drilled 
flow passages. Heating is supplied by means of two electrical band heaters 
clamped to the outside of this cylinder. Heat input is controlled by a commercial 
temperature controller. 

B.2.4 Sample Holder and Diffuser/Thermocouple Disks 

As originally designed the test samples had Viton rubber molded around them. 
Five fine- wire thermocouples were then installed on each face of the regenerator. 
Finally diffuser disks were glued onto both faces, outside the thermocouples. 

This process of vulcanizing Viton around the sample was tried at CSTR, but 
was found to be difficult to predict and control. The process variables (molding 
heat and pressure) had to be modified for each sample. 

Because of this and other problems we devised an alternative process. We 
decided to place the samples in Torlon holders. Torlon (an Amoco Corp. trade- 
name for polyamide- imide) was selected because of its good temperature and 
insulating properties, as well as the availability of a particular grade having ex- 
pansion characteristics almost identical to material of the samples themselves. 
Two separate methods were devised here, one each for the two different types 
of samples which we planned to test. 
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In the case of a sintered rigid material, Torlon holders were heat shrunk 
around the sample. Samples were cooled with liquid nitrogen and the Torlon 
heated in an oven before assembly. 

For stacked screen samples, the screens were first aligned and clamped be- 
tween centers in a lathe. The outside surface of this stack was then coated with 
a high temperature (260 C, 500 F) epoxy, and allowed to cure. Once cured, the 
outside surface of epoxy was machined round. Then several screens were peeled 
off each face of the sample to determine the level of penetration of the epoxy. 
Here instead of a heat shrink process the Torlon holders were fabricated with 
inside diameter O-rings for seeding against the outside of the sample. 

To provide for the thermocouples and the diffuser disks we fabricated reusable 
holders, separate from the test sample. These were fabricated from G7 phenolic. 
Sealing on the face of these next to the heat exchangers was accomplished by 
means of silicone glue. On the faces adjacent to the test sample, sealing was 
provided by O-rings with the grooves cut into the face of the Torlon holders. 
A picture of a diffuser/thermocouple disk and two test samples in their Torlon 
holders is shown in figure B.6. 

B.2.5 Cooler 

The cooler is a shell and tube heat exchanger using water as the coolant. Orig- 
inally type K thermocouples had been selected for measuring the temperature 
rise of the coolant. To improve the measurement of heat rejection we replaced 
these with high accuracy (0.1 C interchangeable) thermistor probes from Omega. 
Here we obtained several probes (for spares) and selected the two best matched 
thermistors from the group. 

To further improve the measurement of heat rejection we replaced the origi- 
nal coolant flow sensor with an Omega model FTB601 having much less pressure 
drop and a lower flow range (linear range 0.1- 2.0 LPM). This flow sensor al- 
lowed us to reduce coolant mass flow significantly and also allowed us to install 
a constant head flow system. 

To solve the problem of inlet coolant (tap water) temperature changes over 
time, we installed a laboratory cooling bath. Cooling water passed through the 
bath in coiled copper tubing, before proceeding to the test rig. 


B.3 Data Acquisition 

The DAS (Data Acquisition System) of the original Pressure Drop Test Rig 
consisted of a Metrabyte analog-to-digital board installed in an IBM compatible 
computer. 

For regenerator testing, we added a submultiplexer board to handle the nu- 
merous additional non-dynamic signals (mostly thermocouples). Here we chose 
(because of cost) a board advertised as equivalent to a Metrabyte submulti- 
plexer, but made by a different manufacturer. 
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Figure B.6: Diffuser /thermocouple disk and test samples in Torlon holders. 
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After much trouble-shooting we tracked erroneous temperature readings to a 
poorly sized resistor in the open thermocouple detect circuits on this board. The 
company acknowledged our finding and sent us replacement resistors. Because 
of time constraints, we did not install these resistors but decided to simply 
disable the detect circuits. 


B.4 Test Procedures 

Pressure-drop testing is relatively straight forward. With only the regenerator 
test sample bolted to the end of the piston cylinder, the rig is charged to pressure 
and the motor is swept over a range of strokes and frequencies, during which a 
number of data points are logged. The important measured quantities are piston 
displacement and cylinder pressure as functions of time. These are integrated 
together to give PV power dissipation, which is the actual basis for friction- 
factor data reduction. During testing, the operator keeps an eye on the piston 
displacement and pressure traces as displayed on an oscilloscope to make sure 
that meaningful data is being logged, which he does by presang the appropriate 
key at the computer console. The process is repeated with at least two working 
gases (helium and nitrogen) to ferret out any gas-property dependence in testing 
or data reduction where none should be. 

Heat-transfer testing is substantially similar. This time, however, the com- 
plete assembly shown in figure B.2 is in place, the heater temperature raised 
and heat rejection to the cooler is added to the important measured quantities. 
The elevated heater temperature poses a tactical difficulty because of time re- 
quired to establish equilibrium temperatures in the entire apparatus. Because 
of this, the rig is generally heated and operated for several hours before log- 
ging data. A change in charge pressure or gas requires additional equilibrating 
time. And during testing the thermal boundary conditions tend to change with 
changing piston stroke and frequency. This requires the operator to log several 
zero-amplitude points over the course of a run so that the background compo- 
nent of heat rejection to the cooler (steady conduction through walls, flanges, 
insulation and the matrix itself) can be calibrated out of the data. 
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C.l Heat-Transfer Model 

The two main mechanisms for axial energy transport in regenerators are en- 
thalpy transport and axial conduction. Therefore, the net cyclic regenerator 
energy flux per unit void area may be expressed as 

qt = qh + qk (C.l) 

where qh is enthalpy transport per unit void area given by 

q h = (cpgT) (C.2) 

and qk is apparent axial conduction per unit void area given by 

% = (C.3) 

The operator { ) indicates time average and k a is the apparent axial conduc- 
tion coefficient (a function of flow variables) comprising molecular conduction, 
matrix solid conduction and thermal diffusion arising from small-scale flow ed- 
dies. The above form for qk is valid since all three effects depend linearly on 
temperature gradient. We assume is roughly constant, and gas and solid 
temperatures are nearly identical so that ^ is representative of the solid tem- 
perature gradient as well. 

The above expression for qh is not yet useful to us because it does not 
explicitly contain the film heat transfer coefficient h for which we are trying to 
solve. However, it is straightforward to make the dependence of qh on h explicit. 
Substituting the identity T ~ (T — T 8 ) + T s into equation (C.2) gives 

qh = (Cp9(r - T,)> + (cpgT e ) (C.4) 

Now, the second time average on the right is negligible for two reasons: (1) the 
variation in T 8 is small and (2) it is 90 degrees out of phase with g — at least for 
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sinusoidal flow. Furthermore, T - T t in the first time average may be replaced 
using the approximate equation 


T.-T = 1 


(C.5) 


\hs J dx 

leaving the following expression for qn 



$ 

ii 

50 Ro 

dT /<A 

dx\h / 

(C.6) 


now, explicitly containing h . 

Equation (C.5) comes from a simplified form of the one-dimensional gas 
energy equation for the porous matrix. The starting point is the complete 
equation in the form 

J^(pe) + ^ ((pe + P)u + q)-Q = 0 ( c - 7 ) 

Matrix-to-gas heat transfer per unit void volume Q can be expressed in terms 
of a film heat transfer coefficient h as 

Q = hs(T, - T) (C.8) 

Ignoring the u 2 / 2 term in e and expressing pe - (c„/ R)P and pe + P = CppT 
gives the alternate form of the energy equation 

W!) ¥ + S ( ‘ wTt ' ,tMr ' TJ "' (C ' 9) 

At this point, some further assumptions simplify the energy equation: 


1. is negligible. 

2. g is uniform in x (although it varies with time). 

&T 

3. The matrix temperature profile is linear, with constant 

4. Apparent conductivity k a is not a function of x. 


Assumption 1 means that the first term in (C.9) is negligible. This is valid 
since there is no significant pressure swing in the proposed rig. Assumption 2 
allows us to factor g from the second term of (C.9). This is valid to the extent 
that assumption 1 holds and temperature at any location in the matrix does 
not fluctuate too much. Assumptions 3 and 4 allow us to neglect for then 


_ 1 . &T . 
With the 


above assumptions, the gas energy equation for the porous matrix 


becomes 


= M T‘-T) 


(C.10) 
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which is easily solved to give equation (C.5) for the gas-to-matrix temperature 
difference. 

Equations (C.3) and (C.6) form the basis for subsequent data modeling. 
They are valid to the extent that the approximations used in obtaining them 
are valid — which some readers are likely to question at this point. While it is 
beyond the scope of the present report, a complete discussion of error estimation 
can be found in the earlier report [3]. 


C.1.1 Dimensionless Forms 


The next step is to re-write the equations for qh and Qk in terms of dimensionless 
groups. Because of the simple expression for qn that will result, we choose to 
represent h in terms of Nusselt number rather than Stanton number. At any 
rate, the two are related by N u = StR^Pr- Starting with equation (C.6) and 
rearranging in terms of Nusselt number (hd^/h), Peclet number (gd^Cp/k) and 
hydraulic diameter (4/s), gives 


9h = 




(C.ll) 


To represent apparent axial conductivity, we choose the dimensionless group- 
ing Nk = k a /k y the ratio of apparent conductivity to gas-molecular conductivity. 
In terms of Nk, (C.3) reduces to 


= (C.12) 

That N u and Nk are good choices for dimensionless groups, is born out by 
the simplicity of the resulting expressions for qh and qk and the common factor 
k occurring in both. 

In our test rig it is the enhancement of axial conduction above molecular 
conduction we are correlating, or qk — qk o- It is not our purpose to measure qk 0 - 
In fact, we have essentially ignored it by calibrating it out of the data acquisition 
process. What we have measured is the total regenerator energy flux less static 
conduction — what we call ql — defined as 


<7t = % + Qk - qko 


(C.13) 


This is because there is no good way to separate qko from the other static 
conduction leaks that must be calibrated out of the test rig. Since, our method 
of measuring regenerator energy flux is not direct but, rather, based on inference 
from cooler heat flux, piston PV power, and various thermal conduction leaks 
through canister and flanges, using qf as the fundamental experimental variable 
just means we lump qko in with the other static conduction leaks. 

Equations (C.ll) and (C.12) enable us to solve for N u and Nk as functions 
of appropriate dimensionless groups of the flow (We will return to this). But 
what are the appropriate dimensionless groups? Dimensional analysis for heat 
transfer problems, in steady and oscillating flow, have led past researchers to 
the following possibilities: 
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• Reynolds number Re 

• Prandtl number P r 

• Valensi number V a 

• Length to diameter ratio Ljdh 

• Mach number 


• Grashof number 


For regenerator flow we categorically reject the Mach and Grashof numbers. 
Mach number is only important when flow velocities approach sonic velocity — 
our flow velocities should never exceed the order of 10 m/s. Grashof number is 
only important when buoyancy forces are significant. 

More tentatively, we reject Ljdh and Valensi number. Ljd h measures the 
ratio of the total matrix flow length to the depth into the matrix required to 
establish velocity and temperature fields. Normally, Ljdh. is sufficiently large 
that we can ignore entrance effects. Valensi number is proportional to the 
square of the ratio of hydraulic diameter to viscous boundary layer thickness in 
oscillating flow. It measures the tendency of fluid velocity near the middle of 
flow channels to lag behind that near the walls in response to oscillating pressure 
gradients. In laminar channel flow, this phenomenon does not seriously affect 
heat transfer until V a approaches about 100. Normally, V a is no more than the 
order of 10 for most regenerators. 

This leaves only Reynolds number and Prandtl number as potential corre- 
lating groups. The question is: Re, P T separately, or together as Peclet number 
p e — R t P r 'f Roughly speaking, R e is inversely related to the thickness of the 
viscous boundary layer while P r is proportional to the ratio of the viscous to 
thermal boundary layer thickness. It is somewhat intuitive, then, to expect film 
heat transfer to correlate with the product of Re and P r . Schlichting ([11], PP- 
298), for example, concludes that for laminar flow over a flat plate 


N u <x 


{ 


JPrSRj 


for P r — 0 
for P r > 0.6 


(C.14) 


The question of correlating in terms of P T is somewhat moot anyway, since past 
evidence suggests that all experiments will be done with gases where P r ~ 0.7. 
At least, if we chose P e to represent our data, then there is some hope that our 
correlations will apply to non P r = 0.7 cases (liquid metals, etc.). 

So we are led to chose correlations for N u and JVj, in terms of Peclet num- 
ber alone, as in correlations (2.6) and (2.7). Of course when extended over a 
range of similar matrix structures, we introduce other relevant dimensionless 
characterizations such as porosity. 
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C.1.2 Separating and qk 

So how is it that equations (C.ll) and (C.12) enable us to solve for N u and 
Nk independently when we are measuring Qh and qk combined together in total 
regenerator energy flux? The completely honest answer is they cannot. If we 
are allowed no further assumptions about qh and qk then there is no way they 
can be separated from q t *. But if we assume only that both N u and Nk scale as 
Pedet number raised to some exponent, and that the exponent is less than one 
and about the same for both, then we can separate them. These assumptions are 
justified by the typical success by others in correlating N u and Nk in exponential 
forms and by arguments for similar scaling of N u and Nk presented in appendix 
F. 

The key observations are the dependence of q h in equation (C.ll) on the 
time-average of P^/N u and of qk in equation (C.12) on the time average of Nk 
only. If N u and Nk both scale with P®, then the qh and qk will scale respectively 
with the time- average of ~ a and P%. So if a < 1, the qk part of total heat flux 
will be relatively more significant at low P c values and the qh part will be more 
significant at high P t values. For example, if a = 1/2 (not too unreasonable) 
then qh will scale as P^ 2 while qk will scale as P^ 2 . It is not to hard to imagine 
then, that testing over a broad enough range of Peclet numbers will allow our 
data reduction process to accurately estimate exponent a. 

C.1.3 Cycle- Averaged Form 

The preceding analysis suggests an alternative method of correlating total re- 
generator energy flux qj. Instead of correlating instantaneous N u and Nk in 
terms of instantaneous Peclet number, it might be sufficient for many purposes 
to correlate ql directly in terms of the peak Peclet value F t m or amplitude, 
presuming it varies sinusoidally over a cycle period. This is plausible because 
we have shown above that ql may be written in the form 

q' t =-k^(f(P e )) (C.15) 

where f(P e ) is some function of Peclet number. If we assume that the time 
average of the function of Peclet number on the right scales with the peak Peclet 
number, then we are led to a correlation of the form (2.4), for the dimensionless 
group 

*< = JW (C- 16 ) 

C.2 Pressure-Drop Model 

Flow friction in the matrix produces pumping loss, which we indirectly measure 
in the form of piston PV power. For pressure drop tests, the sample alone is 
attached to the piston cylinder, without cooler, heater or flow diffusers. Letting 
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w r denote the mean effective pumping power per unit void volume, it follows 
from the usual pumping dissipation arguments that 

-=-(!“) (C - 17) 

where is the regenerator frictional pressure gradient, u is the local void- 
average gas velocity. Expressing the frictional pressure drop in the standard 
Darcy friction factor formulation gives 

dP = (C.18) 

dx dh 2 

where dh is hydraulic diameter and p is density. Substituting this equation into 
the previous one gives 

W r = ^ (fpuM) ( C ' 19 ) 

It is convenient to go one step further and rewrite this in terms of mass flux 
per unit area g = pu, which was one of the key variables used for the previous 
heat-transfer models. 

W r = 

We then correlate / in terms of Reynolds number, as in correlations (2.1) or 

( 2 . 2 ). 

C.3 Model Evaluation 

An important part of the parameter estimation process we are about to describe, 
has to do with actually evaluating the right-hand sides of equations like (C.ll), 
(C.12) and (C.20) for comparison against measured experimental values for the 
left-hand sides — regenerator heat flux and pumping dissipation. This is done 
by numerical simulation of the test rig using a finite difference algorithm based 
on former GLIMPS [5] computational algorithm. Each data point requires a 
separate simulation which is adjusted to match experimental measurables. For 
pressure-drop testing the friction factor in the matrix is adjusted until the piston 
PV power is correct. For heat-transfer testing, both matrix friction factor and 
heat-transfer coefficient are adjusted until PV power and cooler heat rejection 
are correct. The purpose of the simulation is primarily to provide mass flow rate 
through the sample as a function of time. The friction factor and heat-transfer 
coefficients in the simulation have nothing to do with the final correlations 
obtained by fitting the data to equations (C.ll), (C.12) and (C.20). 




Appendix D 

Parameter Estimation 


Two good references on data modeling are: Y. Bard, Nonlinear Parameter 
Estimation [1] and W.H. Press Numerical Recipes [10]. The first book is mostly 
theoretical but very complete. The second book treats data modeling as just 
one topic of the spectrum of numerical arts. It is concise, readable and geared 
toward practicality, with source code listings for all its algorithms. Our data 
reduction mostly followed Numerical Recipes' recommendations. 

In its most abstract form, nonlinear data modeling begins with a model 
equation 

y = f(x, a) (D.l) 

where y represents the dependent experimental variable x represents the in- 
dependent experimental variables and a represents a set of parameters, a* for 
k = 1, . . . , M which enter into the functional expression for / in a nonlinear way 
and are to be determined to best fit a set of data. Equation (D.l) might corre- 
spond to correlations (C.ll), (C.12) or (C.20), where N U1 Nk or / are written 
as expressions involving unknown parameters. 

Given a set of N experiments, indexed by i = 1, . . . , N, one may express the 
measured y* in terms of the model equation (D.l) as 

Vi = /(Xi, a) + Ci (D.2) 

where the e* are the so-called residuals — the differences between the measured 
and theoretical values of y. If the model is good, then the residuals are predomi- 
nantly due to random and normally-distributed measurement errors. Assuming 
that this is the case, and assuming we can estimate in advance the standard de- 
viations <Ji of the residuals, then the maximum likelihood estimate for the model 
parameters a* is found by minimizing the least-squares objective function 

^( a ) = E ( y< ~ - " ~) 2 (D3) 

Even when the residuals are not random, equation (D.3) still makes good intu- 
itive sense as an objective function. 
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The reason for the x 2 notation goes back to probability theory. If the resid- 
uals e, are normally distributed with standard deviation a,, then the minimum 
X 2 (a) (denoted Xmin) turns out to have the well ' known chi-square distribution 
with N-M degrees of freedom. This is in the sense of large numbers of exper- 
imental realizations. We use this observation to check how well our model fits 
the data since the probability that a chi-square distributed variable with N-M 
degrees of freedom exceeds a certain value by chance can be easily computed. 
Thus, given the value Xmi„, one 06111 determine how likely it is. If it is not likely 
at all (probability less than 10" 3 or so), then it must be the case that either (1) 
the model does not fit the data or, (2) the a t were underestimated. Depending 
on degree, alternative (1) is not necessarily bad. For a very careful experiment 
with low measurement error, one expects systematic model errors to dominate. 

To get a more complete picture of the degree to which the model fits the data, 
we plot the normalized residuals e./a, vs. peak Reynolds number. Theoretically 
for a good data fit, this plot should be randomly distributed with standard 
deviation 1. Random fluctuations larger than this signify underestimated cr t ; 
non-random trends signify an imperfect model. 

Once we have determined the maximum likelihood parameters a k and con- 
cluded that they do indeed fit the data, we have yet another question we must 
' ask: To what degree of accuracy have we determined the a*. In other words: 
Might some other nearby values of the a k fit the data nearly as well? Or, what 
is the sensitivity of the model to changes in the a* parameters? 

The following ideas are taken pretty much from section 14.5 of Numerical 
Recipes. Essentially we are trying to determine analytically the spread of model 
parameters a that would result from a large number of experimental realizations 
of our data set. We are focusing on just one component a k at a time and want 
to predict its distribution in the universe of experimental realizations. 

First, what is the general relationship between changes in x 2 and changes 
in a? Evidently, this depends on the functional form of the model equation 
(D.l). Expanding the x 2 objective function in a Taylor series at the minimizer 
ao, the first-derivative terms drop out, because the gradient is zero, leaving the 
variation in x 2 dependent on the second- derivative terms. This leaves, for any 
small change fia in the model parameters, 

Ax 2 « ^a -[a] -6a (D-4) 


where matrix [a], defined by 

_ 1 ( d/ ( x t ? a ) df ( x * 2 a ) ^ (D.5) 

afcZ ^ cr ? \ da k dai ) 

is an approximation to 1/2 the Hessian (second derivative matrix) of the x 2 
objective function. As it turns out, and this will be important later, [a] 
is the so-called covariance matrix [C], calculated in the Levenberg-Marquardt 

Algorithm. , 

Now we restrict ourselves to a fixed k and a curve C m M-dimensional 

parameter space defined as the locus of the following points: 
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• a,k a random variable from the universe of experimental realizations, and 

• the other parameters (o/ for l ^ k) adjusted (as a function of a *) to 
minimize x 2 > which is evaluated using a fixed data set. 

It turns out that the distribution of x 2 for the points on C is chi-square with 
one degree of freedom. 

Along C y by definition, all but the fc-th component of the gradient of x 2 is 
zero. But the gradient of x 2 may be written 

Vx2 ^^£' 5 ° i = 2[a] ' 5a (D6) 

so that for 5a tangent to C at the minimizer ao 

[a] * 5a = bS k (D.7) 

where 5* is the vector with one in its k - th component and zero elsewhere and b 
is some constant. Substituting (D.7) into (D.4) gives 

Ax 2 » 5a • b5jfe = b6a * (D.8) 

along C. But also, (D.7) may be solved to give 

5a = [<x]~ l * &5fc = [C] * b6k (D.9) 

or, in particular, 

6ak = bCkk (D.10) 

Eliminating b between equations (D.8) and (D.10) gives 

6a k = y/A (D.ll) 

which holds so long as 5a is tangent to C at ao. 

Since we know the distribution x 2 along C is chi-square with one degree of 
freedom we are in business. The p-confidence interval for parameter a* may be 
expressed __ 

6<Lk = ±y/A py/Ckk (D.12) 

where A p , is the number for which the probability is p that a chi-square variable 
with one degree of freedom is less than A p . A p can be found in tables for several 
values of p, or can be computed numerically as a function of p. For example, 
with 68.3 percent confidence (p — 0.683), A p — 1, with 90 percent confidence 
A p = 2.71, etc. 

Equation (D.12) is readily evaluated, but we must keep in mind that its 
validity depends on a good fit of the model to the data. If the value Xmin f Qun d 
in the function minimization process is too large (unlikely) then our assumption 
that it is distributed as a chi-square variable with N — Af degrees of freedom 
is probably wrong. In this case our A p estimate (made as above) would be 
underestimated and our confidence intervals would be too small. 
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Software 


Data reduction is a multi-step process where hardware data is converted into 
final heat-transfer or pressure-drop correlations. Along the way, a number of 
computer files come into being. These start out with run-date encoded names 
and progress to more descriptive names later on. These names, and other im- 
portant information, are hand recorded on log sheets for future reference. 

The individual computer programs that make data reduction possible are: 


Heat-Transfer Reduction 


program 

RRLSEL 

input 

interactive 

ouipuu 

console 

urau 

select sample length 

REGENR1G 

RR file 

.RAW file 

source data 

RRHSCRN 

.RAW file 

.SCN file 

data screening 

RRHTRANS 

.SCN file 

.DRV file 

data processing and translation, including 
numerical modeling 

COMBINE 

.DRV rile 

.DRV file 

combines files together 

RRVIEW 

.DRV file 

.DRV file 

data viewing, sorting, selecting 

RRHFUDGE 

.DRV file 

.DRV file 

data correction 

RRH1, ... 

.DRV file 

console 

data modeling 



Pressure-Drop Reduction 

program 

input 

output 

description 

REGENRIG 

.RR file 

.RAW file 

source data 

RRPSCRN 

RRPTRANS 

.RAW file 
.SCN file 

.SCN file 
.DRV file 

data screening 

data processing and translation, including 



numerical modeling 

COMBINE 

.DRV rile 

.DRV file 

combines files together 

RRVIEW 

.DRV file 

.DRV file 

data viewing, sorting, selecting 

RRPFUDGE 

.DRV file 

.DRV file 

data correction 

RRP1, ... 

.DRV file 

console 

data modeling 
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Except for COMBINE, which is a DOS batch file, these programs are source- 
coded in Microsoft Pascal and run under MS-DOS. 

The following account is a rough outline of the data reduction process. In- 
dividual programs, previous report [3] and a number of unpublished personal 
memoranda contain more detailed comments. 

Size Sample Length Prior to testing, program RRLSEL estimates the best 
sample length for accurate data reduction over the broadest range of Reynolds 
number. 

Log Data During testing, program REGEN RIG logs a number of data points 
producing a .RR data file in binary format. 

Convert to Raw Format Later, program REGENRIG reads the .RR file as 
input and writes a .RAW data file by the same name as output. The .RAW file 
has ASCII format and can be modified with a text editor to delete bad data 
points by hand, etc. 

- Screen Data Program RRPSCRN or RRHSCRN screen .RAW files, removing 
bad data points and performing elementary processing (such as determining 
static conduction from zero-amplitude points), producing .SCN files. Often 
run in batch mode on multiple files using RUNSP.BAT or RUNSH.BAT. 

Simulation Program RRPTRANS or RRHTRANS read .SCN files as input, call 
on numerical simulation and optimization modules to duplicate the test point, 
producing binary .DRV files. Often run in batch mode on multiple files using 
RUNTP.BAT or RUNTH.BAT. 

Combine Batch file COMBINE.BAT concatenates individual .DRV files into a 
single .DRV file, usually to combine runs with the same test sample taken on 
different days. Done with the DOS COPY /B command, or similar. 

Select Program RRVIEW reads .DRV files enabling the user to individually 
view, discard and sort data points, if necessary. 

Correct Programs RRPFUDGE and RRH FUDGE allow the user to modify in- 
dividual fields for each data point of a .DRV file, if necessary. Often run in batch 
mode on multiple files using RUNFP.BAT or RUNFH.BAT. 

Data Modeling Programs RRP1, RRH1, ..., read .DRV files, fitting coef- 
ficients of correlating expressions to data, producing tabular output (in file 
RRMODEL.OUT) and graphical output. Graphical output can be printed with 
a DOS screen printing utility. 
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N u and Nj- Physics 


In correlations (2.6) and (2.7) we assumed that Nusselt number N u and ax- 
ial conductivity enhancement ratio N k both scaled in proportion to Reynolds 
number raised to the same power. To justify this in porous flows, we can make 
use of the same methods developed for study of turbulent flows — namely, 
that enhancements over molecular thermal diffusion are understood in terms 
of fluctuating parts of velocity and temperature fields. In two dimensions, for 
example, u = u + u' and v = v + v' are the x and y velocity components and 
T = T + T' is temperature, the — operator denoting the local steady part and 
' denoting the fluctuating part. Even though flow through porous materials 
is not strictly turbulent, it does have small fluctuating components superim- 
posed on the steady fluid variables and therefore can be understood in the same 
terms. Irregularities present in porous flow include flow separations caused by 
abrupt expansions, recirculations, vortex shedding, channeling and even true 

turbulence at high enough Reynolds numbers. 

In our case of oscillating flow, the steady part of a fluid variable refers to 
its ensemble average in the sense of a large number of instances at the same 
position and crank angle. The fluctuating part is the instantaneous deviation 
from the steady part. But even the steady part of a flow variables is on too 
fine a scale to be resolved in a typical computational solution. So we will talk 
of yet another average value, the mean-flow value, which is the steady part 
spatially- averaged over a region comparable in size to the computational grid. 
Steady and fluctuating parts are denoted with overbars and primes as above. 
We have no special notation for mean-flow values, preferring to just state the 
fact in English when required. 

The key phenomenon is the diffusion of thermal energy brought about by 
fluctuating velocity components in the presence of a temperature gradient. This 
occurs at a scale not resolved by the usual computational grid, hence the need 
for empirical notions such as enhanced axial conductivity and film heat transfer 
coefficient, non- dimensionally represented in terms of axial conductivity ratio 

JV fc and Nusselt number N u . . 

Chapter XXIII of Schlichting [11] shows that, for an incompressible fluid, 
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the enhanced heat flux per unit area passing through planes normal to the x 
and y directions is 

qi = pCpU'T' (F.l) 

<h = pep^T (F.2) 

It is easy to see that these are just those parts of standard enthalpy-flux inte- 
grals not resolved by the steady flow variables. Less clear is why velocity and 
temperature fluctuations should be correlated in the presence of a temperature 
gradient. To understand why, assume there is a steady temperature gradient, 
say 7 ^ 0 , and let your point of observation be centered on a fixed-in-space 
microscopic volume element on the same order of size as the smallest flow eddies. 
Into this volume element comes a fluid eddy perturbed from the steady velocity 
by amount u '. At this point it helps to think of u ' as just one component of the 
full spectrum of velocity fluctuation, say v! = A^cosut. Because of the temper- 
ature gradient, the eddy temperature differs from the steady temperature of the 
surroundings by an amount |^Ax, where Ax = f t u f — (A^/uS) sin wt is a sort 
of mixing length representing the present displacement of our eddy compared 
to the steady flow. This temperature fluctuation does not contribute to en- 
hanced heat flux because Ax and v! are 90 degrees out of phase, canceling over 
time. However, it does drive a flow of heat, by molecular conduction, to or from 
the local surroundings, roughly Q a f-Ax, in turn producing an additional 

temperature fluctuation J t Q oc || J t Ax oc ~§|(A^/^ 2 )cosu;*. Evidently this 
temperature fluctuation does correlate with u f and produces the enhanced heat 
flux. It is seen that enhanced heat flux is similar to a shuttle heat transfer mech- 
anism operating between each flow eddy and its ensemble- average surroundings. 

At this point we conclude that the enhanced heat flux vector at any point 
may be represented in the form q = k t VT , where k t , the eddy conductivity, 
depends only on the nature of the local velocity fluctuations and VT is the 
local steady temperature gradient. Generally we should expect that k t varies 
from point to point and time to time in the porous flow Held, and also varies 
with the direction of VT at any given point. But we will get nowhere thinking 
like this so we will now make an assumption about k t which will turn out to 
be important: we assume that kt is a function only of the mean-flow Reynolds 
number and that it varies with Re in the same way at each point of the porous 
flow field. This would be true, for example, if k t were everywhere proportional 
to R% for some exponent a. 

In order to talk next about enhanced axial conduction and film heat transfer, 
we need to introduce some notation regarding the regions within which these 
phenomena are taking place. In both cases we will be referring to a volume 
element C of a porous matrix — a cube aligned with the coordinate axes whose 
side dimension dx is much larger than the pore size of the matrix but still 
small compared to the computational grid, so that mean-flow variables may be 
considered uniform within C, The void volume within C is denoted V, with 
volume V . And the surface of V, which is mainly the surface of the gas-solid 
interface within C, is denoted 5, with area S'. 
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We will also require another bold assumption: that the steady fluid tem- 
perature gradient may be decomposed everywhere into a part that drives en- 
hanced axial conduction and a part that drives film heat transfer. In notation 
VT = VT k + VT*. Another way of looking at this is to suppose that the 
enhanced heat flux vector at every point of the flow field contributes partly to 
enhanced axial conduction and partly to film heat transfer. We just imagine 
that said decomposition exists implicitly, and is uniquely determined, not that 
we can expressly evaluate it. 

F.l Axial Conductivity Ratio 

To best understand enhanced axial conduction, imagine our volume element C 
with fixed temperature boundary conditions imposed on the axially opposite 
faces so that there is a mean-flow axial temperature gradient set up. This tem- 
perature gradient will impress itself on the fluid everywhere within void volume 
V, with local steady variations to some degree. These local variations may be 
supposed significant in some cases, but VT* should always scale in proportion to 
the mean-flow temperature gradient. Clearly, the previously described enhanced 
heat flux mechanism will be active within void volume V and the individual en- 
hanced heat flux vectors will sum together into an effective mean-flow enhanced 
axial heat flux for the entire volume. It is this mean-flow enhanced axial heat 
flux, divided by that expected from molecular conduction alone, which is our 
dimensionless group N k . Evidently, N k depends only on some weighted average 
of kt over V and with our previous assumption as to the nature of kt, we con- 
clude that N k must vary with mean-flow Reynolds number in the same way as 
k t does. 

F.2 Nusselt Number 

Relating Nusselt number to k t is also possible. Consider again our volume 
element C, this time focusing on film heat transfer. Taking first the usual 
macroscopic engineering point of view, the heat flow between the fluid and 
solid is formulated as 

Q = hS&Tf (F-3) 

where h is a film heat transfer coefficient and AT/ is the mean-flow to solid 
surface temperature difference. (Here, we have ignored that relatively-small 
part of S formed by the outside boundary of V that falls within the fluid.) But 
this same heat flow is also the surface integral over S of the normal heat flux, 
or f 

Q — J q ' (F-4) 

where n is the unit surface normal vector. To get away from the surface where 
molecular conductivity dominates we need to translate this into a volume inte- 
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gral using the divergence theorem which says, in this case: 

Q= [ V-qdv (F.5) 

Jv 

Presuming eddy conduction dominates (otherwise we would be in the limiting 
range of constant N u — Auo)> and using our above formulation for q, we can 
now write 


Q= f V- k t VT h dv= f (k t V 2 T h + VJ* • VT fc ) dv (F.6) 

Jv Jv 

Equating the right sides of equations (F.3) and (F.6) allows us to solve for film 
coefficient ft, and then Nusselt number as 

= + (F.7) 

where is hydraulic diameter and k is molecular conductivity. If we now make 
the additional assumption that VT& (and therefore V 2 7\) scales everywhere 
in proportion to the film temperature difference AT/, we conclude that Nusselt 
number is proportional to a weighted average over V of k t and Vfc t . If k t varies 
only with R<>, then we might as well assume that Vfct does too, so we conclude 
that Nk must vary with mean-flow Reynolds number in the same way as k t does. 


F.3 Similarity 

Based on the above discussion we find it arguable that both axial conductivity 
ratio Nk and Nusselt number N u vary with mean-flow Reynolds number in the 
same way as does eddy conductivity k t . The main simplifying assumptions 
leading to this conclusion were that: 

• Eddy conductivity k t is the principal phenomenon responsible for en- 
hanced mean-flow axial conduction and film heat transfer, or at least for 
their variations with mean-flow Reynolds number. 

• This eddy conductivity varies with mean-flow Reynolds number, the same 
way at all points of the flow field and in all directions. 

• Local steady fluid temperature gradient may be decomposed into an en- 
hanced axial conduction part and a film heat transfer part, the former 
scaling in proportion to the mean-flow axial temperature gradient and the 
latter in proportion to the film temperature difference. 

Therefore, to the extent this is all true, we are justified in Simultaneously re- 
ducing Nk and N u in the forms of correlations (2.6) and 2.7). Further evidence 
that this is a reasonable thing to do comes in the form of the relatively tight 
confidence intervals obtained for the resulting parameter estimates. 
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